#include <stdio.h>
#include <stdlib.h>
#include <math.h>
//#include <atlas/clapack.h>
//#include "clapack.h"
#include "FDmatrix.h"
#include "materialfunction.h"


void createFDmatrix(int iter_z,double dt,double dz, double material[][3],
                    double phi[],double k0[],double Hpicar[],double H0[],double rhomx[][3], double A[], double b[],int lda){


    double K1,K2,K3,Kmid,Kup,Klow,K1up,K2up,K3up,K1low,K2low,K3low,Kfor,Kback;
    double c1,c2,c3,c;
    double theta1,theta2,theta3,theta;
    double theta01,theta02,theta03,theta0;
    double PHI;
    int i;
    for(i = 0;i<iter_z;i++){
        if (i == 0){
            A[0] = 0;
            // dirichlet boundary

            A[1] = 0;
            A[2] = 1;
            A[5] = 0;
            b[0] = 0;
        }else if( i > 0 && i < (iter_z-1)){

            K1       = permeability(material,k0,Hpicar[i],1);
            K2       = permeability(material,k0,Hpicar[i],2);
            K3       = permeability(material,k0,Hpicar[i],3);
            Kmid     = rhomx[i][0]*K1+rhomx[i][1]*K2+rhomx[i][2]*K3;
            /*printf("K1 %e \n",K1);
            printf("K2 %e \n",K2);
            printf("K3 %e \n",K3);
            printf("Kmid %e \n",Kmid);*/

            K1up     = permeability(material,k0,Hpicar[i+1],1);
            K2up     = permeability(material,k0,Hpicar[i+1],2);
            K3up     = permeability(material,k0,Hpicar[i+1],3);
            Kup      = rhomx[i+1][0]*K1up+rhomx[i+1][1]*K2up+rhomx[i+1][2]*K3up;
            /*printf("K1up %e \n",K1up);
            printf("K2up %e \n",K2up);
            printf("K3up %e \n",K3up);
            printf("Kup %e \n",Kup);*/

            K1low    = permeability(material,k0,Hpicar[i-1],1);
            K2low    = permeability(material,k0,Hpicar[i-1],2);
            K3low    = permeability(material,k0,Hpicar[i-1],3);
            Klow     = rhomx[i-1][0]*K1low+rhomx[i-1][1]*K2low+rhomx[i-1][2]*K3low;
            /*printf("K1 %e \n",K1low);
            printf("K2 %e \n",K2low);
            printf("K3 %e \n",K3low);
            printf("Klow %e \n",Klow);*/

            Kfor     = 0.5*(Kmid+Kup);
            Kback    = 0.5*(Kmid+Klow);

            theta1   = moisturecontent(material,Hpicar[i],1);
            theta2   = moisturecontent(material,Hpicar[i],2);
            theta3   = moisturecontent(material,Hpicar[i],3);
            theta    = rhomx[i][0]*theta1+rhomx[i][1]*theta2+rhomx[i][2]*theta3;
            /*printf("th1 %e \n",theta1);
            printf("th2 %e \n",theta2);
            printf("th3 %e \n",theta3);
            printf("th %e \n",theta);*/


            theta01  = moisturecontent(material,H0[i],1);
            theta02  = moisturecontent(material,H0[i],2);
            theta03  = moisturecontent(material,H0[i],3);
            theta0   = rhomx[i][0]*theta01+rhomx[i][1]*theta02+rhomx[i][2]*theta03;
            /*printf("th01 %e \n",theta01);
            printf("th02 %e \n",theta02);
            printf("th03 %e \n",theta03);
            printf("th0 %e \n",theta0);*/

            c1       = conductivity(material,Hpicar[i],1);
            c2       = conductivity(material,Hpicar[i],2);
            c3       = conductivity(material,Hpicar[i],3);
            c        = rhomx[i][0]*c1+rhomx[i][1]*c2+rhomx[i][2]*c3;

            PHI      = rhomx[i][0]*phi[0]+rhomx[i][1]*phi[1]+rhomx[i][2]*phi[2];           /*printf("c1 %e \n",c1);
            printf("c2 %e \n",c2);
            printf("c3 %e \n",c3);
            printf("c0 %e \n",c);*/



            A[lda*i]     =  0;
            A[(lda*i)-1] = -Kback / pow(dz,2);                    //A(i,i-1)   //A(i,i)
            A[(lda*i)+2] =  PHI*c/dt + (Kfor+Kback) / pow(dz,2);   //A(i,i+1)
            A[(lda*i)+5] = -Kfor /  pow(dz,2);
            b[i]         =  PHI*Hpicar[i]*c/dt + (Kfor-Kback)/dz + PHI*(theta0 - theta)/dt;
        }else{//i <iter_z

            //Neumann boundary
            K1up     = permeability(material,k0,Hpicar[i-1]-2*dz,1);
            K2up     = permeability(material,k0,Hpicar[i-1]-2*dz,2);
            K3up     = permeability(material,k0,Hpicar[i-1]-2*dz,3);
            Kup      = rhomx[i-1][0]*K1up+rhomx[i-1][1]*K2up+rhomx[i-1][2]*K3up;

            K1       = permeability(material,k0,Hpicar[i],1);
            K2       = permeability(material,k0,Hpicar[i],2);
            K3       = permeability(material,k0,Hpicar[i],3);
            Kmid     = rhomx[i][0]*K1+rhomx[i][1]*K2+rhomx[i][2]*K3;

            K1low    = permeability(material,k0,Hpicar[i-1],1);
            K2low    = permeability(material,k0,Hpicar[i-1],2);
            K3low    = permeability(material,k0,Hpicar[i-1],3);
            Klow     = rhomx[i-1][0]*K1low+rhomx[i-1][1]*K2low+rhomx[i-1][2]*K3low;

            Kback    = 0.5*(Kmid+Klow);
            Kfor     = 0.5*(Kmid+Kup);

            theta1   = moisturecontent(material,Hpicar[i],1);
            theta2   = moisturecontent(material,Hpicar[i],2);
            theta3   = moisturecontent(material,Hpicar[i],3);
            theta    = rhomx[i][0]*theta1+rhomx[i][1]*theta2+rhomx[i][2]*theta3;

            theta01  = moisturecontent(material,H0[i],1);
            theta02  = moisturecontent(material,H0[i],2);
            theta03  = moisturecontent(material,H0[i],3);
            theta0   = rhomx[i][0]*theta01+rhomx[i][1]*theta02+rhomx[i][2]*theta03;

            c1       = conductivity(material,Hpicar[i],1);
            c2       = conductivity(material,Hpicar[i],2);
            c3       = conductivity(material,Hpicar[i],3);
            c        = rhomx[i][0]*c1+rhomx[i][1]*c2+rhomx[i][2]*c3;

            PHI      = rhomx[i][0]*phi[0]+rhomx[i][1]*phi[1]+rhomx[i][2]*phi[2];


            A[lda*i]     =  0;
            A[(lda*i)-1] = -Kback / pow(dz,2)-Kfor / pow(dz,2);                    //A(i,i-1)
            A[(lda*i)+2] =  PHI*c/dt + 2*Kback / pow(dz,2);      //A(i,i)
            A[(lda*i)+3] =  0;                  //A(i,i+1)

            b[i]         = PHI*Hpicar[i]*c/dt + PHI*(theta0 - theta)/dt+(Kfor-Kback)/dz+2*Kfor/dz;
        }


    }

}
